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ABSTRACT 

We study transients produced by equatorial disk-like outflows from catastrophically 
mass-losing binary stars with an asymptotic velocity and energy deposition rate near 
the inner edge which are proportional to the binary escape velocity v eS c- As a test case, 
we present the first smoothed-particle radiation-hydrodynamics calculations of the 
mass loss from the outer Lagrange point with realistic equation of state and opacities. 
The resulting spiral stream becomes unbound for binary mass ratios 0.06 < q < 0.8. For 
synchronous binaries with non-degenerate components, the spiral-stream arms merge 
at a radius of ~ 10n, where a is the binary semi-major axis, and the accompanying 
shock thermalizes about 10% of the kinetic power of the outflow. The mass-losing bi¬ 
nary outflows produce luminosities reaching up to ~ 10 6 L Q and effective temperatures 
spanning 500 < T e g < 6000 K, which is compatible with many of the class of recently- 
discovered red transients such as V838 Mon and V1309 Sco. Dust readily forms in 
the outflow, potentially in a catastrophic global cooling transition. The appearance of 
the transient is viewing angle-dependent due to vastly different optical depths parallel 
and perpendicular to the binary plane. We predict a correlation between the peak 
luminosity and the outflow velocity, which is roughly obeyed by the known red tran¬ 
sients. Outflows from mass-losing binaries can produce luminous (10 s L 0 ) and cool 
(T e ff < 1500 K) transients lasting a year or longer, as has potentially been detected by 
Spitzer surveys of nearby galaxies. 

Key words: Binaries: close — binaries: general - stars: evolution — stars: winds, 
outflows 


1 INTRODUCTION 

A fraction of binary stars evolve through at least one phase 
of dramatic angular momentum loss, accompanied by a dras¬ 
tic reduction of the semi-major axis and envelope mass ejec¬ 
tion. The generally accepted framework for this stage of bi¬ 
nary evolution is the “common envelope”, when the drag 
forces on the binary orbiting in an non-corotating envelope 
lead to tightening of the orbit and ejection of the envelope 
(e.g. Paczynski 1976; Iben & Livio 1993; Taam & Sandquist 
2000; Ivanova et al. 2013a). In some cases the binary may 
survive this phase with modified orbital parameters, be¬ 
coming visible as a cataclysmic variable, X-ray binary 
or recycled pulsar (e.g. Meyer & Meyer-Hofmeister 1979; 
Bhattacharya & van den Heuvel 1991). However, in other 
cases the two stars may merge completely into a single ob¬ 
ject. Perhaps the best known example is the progenitor of 
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SN1987A, which is surrounded by a triple-ring structure, 
presumably ejected during such a merger (Chevalier & Soker 
1989; Hillebrandt & Meyer 1989; Podsiadlowski et al. 1990; 
Podsiadlowski 1992; Morris & Podsiadlowski 2007, 2009). 
Alternatively, the tightened orbit might lead to a delayed 
merger driven by gravitational wave emission or by an action 
of distant companions through the Kozai-Lidov process (e.g. 
Thompson 2011; Katz & Dong 2012; Pejcha et al. 2013), ul¬ 
timately involving the extreme physics of compact objects. 
This channel includes Type la supernovae important for 
distance measurements in cosmology (e.g. Iben & Tutukov 
1984; Webbink 1984; Weinberg et al. 2013), and neutron 
star and black hole mergers, which are promising sources 
for detecting gravitational waves (e.g. Abadie et al. 2010; 
Metzger et al. 2010). 

Despite its importance to many areas of astro¬ 
physics, observational tests of the common envelope frame¬ 
work have primarily been indirect and based on compar¬ 
ing observed post-common envelope binaries to theoret- 
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ical predictions, often aiming to constrain the parame¬ 
terized common envelope efficiency (Webbink 1984). Re¬ 
cently, this dynamic and rapid phase of binary evolu¬ 
tion has been associated with transients that have lu¬ 
minosities in between classical novae and supernovae, 
and red color (e.g. Munari et al. 2002; Bond et al. 
2003; Soker & Tylenda 2003, 2006; Tylenda & Soker 2006; 
Kulkarni et al. 2007; Thompson et al. 2009; Tylenda et al. 
2011, 2013; Ivanova et al. 2013b). These transients, which 
include V838 Mon, M85 OT, M31 RV, and V4332 Sgr, 
are commonly named ‘luminous red novae’, ‘interme¬ 
diate luminosity optical transients’, or ‘mergebursts’ 
(Soker & Tylenda 2006). Hereafter, we describe all of these 
events using the generic term ‘red transients’ (RT), because 
some members of this class are as luminous as, or even dim¬ 
mer than, classical novae. The association between RT and 
binary stars has been strengthened by the pre-outburst pho¬ 
tometry of V1309 Sco, which reveals an eclipsing contact bi¬ 
nary with a rapidly decreasing orbital period (Tylenda et al. 
2011). It has been proposed that RT optical emission re¬ 
sults from the release of thermal and hydrogen recombina¬ 
tion energy from a shell heated and launched by the binary 
interaction, similar to a shock passing through the hydro¬ 
gen envelope of a Type II-Plateau supernova and leading to 
relatively cool effective temperatures (Ivanova et al. 2013b). 
RT occur every ~ 2 years in Milky Way-like galaxies, a rate 
which is broadly consistent with binary population synthesis 
models (Kochanek et al. 2014). 

RT exhibit considerable diversity, with peak luminosi¬ 
ties ranging from 1O 3 -1O 4 L 0 for V4332 Sgr and OGLE 2002- 
BLG-360 (Tylenda et al. 2013) to ~ 1O 6 L 0 for V838 Mon, 
or up to 10^ L 0 if objects such as SN2008S and NGC300 
OT are included. Explaining the latter two events also re¬ 
quires an unobserved shock breakout with a luminosity of 
~ 1O 1O L 0 (Kochanek 2011a). Note also that the luminos¬ 
ity of the “Rosetta stone” event V1309 Sco was a factor of 
~ 10 lower than V838 Mon and similar objects, and factor of 
~ 100 lower than SN2008S and NGC300 OT. The outburst 
durations of RTs range from about 20 days for V4332 Sgr 
(Martini et al. 1999) to roughly a year or longer for CK Vul 
and OGLE 2002-BLG-360 (Kato 2003; Tylenda et al. 2013). 

Fundamental differences also exist among the RT pro¬ 
genitors, when such identifications are possible. The pro¬ 
genitor of V1309 Sco was a binary containing an K-type 
primary (Tylenda et al. 2011), M31 LRN 2015 harbored a 
2-4 Mq evolved star (Dong et al. 2015), V838 Mon is known 
to have B3 main-sequence progenitor (Munari et al. 2005; 
Tylenda et al. 2005b), while the progenitor of SN2008S 
was a dust-obscured extreme asymptotic giant branch star 
(Prieto et al. 2008; Thompson et al. 2009). The progenitor 
mass and luminosity appears to be correlated with the peak 
luminosity of the transient (Kochanek et al. 2014). 

The light curves also display dramatic differences: some 
events show relatively smooth initial declines (SN2008S, 
NGC300 OT, V4332 Sgr), while others show multiple peaks 
(V838 Mon, OGLE 2002-BLG-360, CK Vul), or evidence for 
sudden dust formation (M31 LRN 2015, SN2008S, NGC300 
OT). V1309 Sco exhibited a slow rise to maximum, which 
lasted ~ 150 days - a relatively unique feature among the 
RT, shared only by OGLE 2002-BLG-360. However, a simi¬ 
lar feature might have been missed in other objects. Finally, 


the RT also display distinct differences in their spectroscopic 
evolution. 

The rich phenomenology of RT warrants a detailed anal¬ 
ysis of the observable signatures of the dynamical phases 
of binary evolution. By investigating the signatures of each 
process contributing to this evolution individually, we may 
hope to identify those discriminating characteristics which 
ultimately allow the observed diversity of RT to be mapped 
into individual astrophysical processes. More ambitiously, by 
linking the remnants of the observed transients with known 
classes of objects, RT could one day serve as probes of the 
physics of the progenitors and their evolutionary pathways. 
In this paper, we present the first steps in this direction. 

One important channel to remove angular momentum 
from a binary star and shrink its orbit is through mass loss 
from the outer Lagrangian (L 2 ) point, a process which is effi¬ 
cient due to the long lever arm of the L 2 point (e.g. Pribulla 
1998). For L 2 mass loss to occur, both stars need to fill their 
Roche lobes and their common photosphere must expand 
sufficiently. This condition can be achieved by an evolution¬ 
ary change in the structure of one of the stars, which initiates 
a rapid, dynamical mass transfer (Webbink 1976). The mass 
transfer rate will accelerate if there is positive feedback due 
to the adiabatic response of a donor if the latter possesses 
a deep convective envelope (Paczynski & Sienkiewicz 1972; 
Hjellming & Webbink 1987; Soberman et al. 1997; Ge et al. 
2010, 2015; Woods & Ivanova 2011; Pavlovskii & Ivanova 
2015). The secondary component cannot thermally accom¬ 
modate the incoming material, swells up and fills its own 
Roche lobe. Once the common surface reaches L 2 as the 
two stars get closer by mass transfer, mass and angu¬ 
lar momentum loss from L 2 will commence (e.g. Webbink 
1977). Alternatively, the binary can be brought into con¬ 
tact via a tidal instability (e.g. Rasio & Shapiro 1995; 
Eggleton & Kiseleva-Eggleton 2001), which also eventually 
results in dynamically unstable mass transfer and L 2 mass 
loss (Tylenda et al. 2011; St§pien 2011). 

Matter leaves the L 2 point in a spiral stream, which is 
torqued and accelerated by the time-changing gravitational 
force of the central binary. Shu et al. (1979) shows that for 
binary mass ratios 0.064 < q < 0.78 the stream becomes un¬ 
bound and escapes from the system. For q outside of this 
range, the stream remains bound and instead accumulates 
in a mass-loss ‘ring’, which continues to interact with the 
binary. The rest of this paper will focus predominantly on 
the unbound L 2 outflows, deferring a discussion of the ob¬ 
servable signatures of the mass-loss ring to future work. 

Spiral streams leaving the binary spread radially with 
time, causing them to eventually collide and dissipate a frac¬ 
tion of their kinetic energy into heat and radiation by shocks. 
Although the merging of the spiral arms has long been ex¬ 
pected (Shu et al. 1979), including the potential for shocks 
(Leibowitz & Mazeh 1981), the importance of L 2 mass loss 
events as transient electromagnetic sources has thus far been 
neglected. In this paper, we show that the L 2 binary outflows 
can shine with high luminosities and low effective tempera¬ 
tures, similar to the observed properties of the RTs. 

L 2 mass loss may be of particular relevance for inter¬ 
preting the light curve of V1309 Sco, which exhibits three 
unexpected features, as pointed out in Pejcha (2014). First, 
the decrease in the orbital period of the contact binary di¬ 
rectly observed in the ~ 2000-day light curve prior to the 
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main outburst is rapidly accelerating; the second and higher 
period derivatives are incommensurably higher than the de¬ 
cay timescale P/P. Second, the observed gradual changes 
in the phased binary light curve can only be understood if 
the binary is obscured anisotropically from one direction. 
Finally, after the contact binary ceased to be visible, its 
gradual brightening to maximum lasted ~ 150 days - much 
longer than the dynamical timescale of the binary given by 
its orbital period of P » 1.4 days. Pejcha (2014) argued that 
these three features are explained by L 2 mass loss, which 
naturally ties together the period change, the obscuration 
of the binary star, and the brightening due to expansion of 
the pseudophotosphere in the outflow as the mass-loss rate 
runs away. Although we defer a more thorough argument 
to future work, V1309 Sco strongly indicates that at least 
in some systems L 2 mass loss runs away rather slowly, over 
many hundreds or thousands of binary dynamical timescales. 
This justifies approximating the mass loss by a sequence of 
steady-state models. 


Our results indicate that L 2 mass loss leads to a 
universal configuration defined by a disk-like equatorial 
outflow, with energy deposited by colliding spiral streams 
at the inner edge. Both the energy deposition rate and the 
asymptotic velocity are proportional to the instantenous 
escape velocity from the binary. Similar outflow configu¬ 
rations are expected to occur also during some common 
envelope scenarios, when the released binary orbital energy 
leads to envelope heating and ejection (e.g. Paczynski 1976; 
Webbink 1984; Iben & Livio 1993; Taam & Bodenheimer 
1989; Livio et al. 1990; Sandquist et al. 1998; 
Taam & Sandquist 2000; Morris & Podsiadlowski 2006; 
Ricker & Taam 2012; Passy et al. 2012; Ivanova et al. 
2013a; MacLeod & Ramirez-Ruiz 2015). Although the 
details of this process are not fully understood, the energy 
deposition rate and asymptotic ejection velocity likely scale 
with the orbital velocity of the binary. Our calculations are 
thus relevant for predicting the observable signatures of this 
considerably wider class of models, and we use L 2 mass loss 
as an example of a relatively clean link between the binary 
properties and the energy deposition rate. 


In this paper, we present the first radiation- 
hydrodynamics calculations of L 2 mass loss, focusing on the 
observational appearance of these events. In Section 2, we 
describe the numerical setup of our calculations. In Sec¬ 
tion 3, we review the semi-analytic theory of the L 2 spiral- 
stream outflow and analyze the numerical results. We also 
constrain the transport of the angular momentum from the 
binary to the L 2 outflow, and estimate its gravitational and 
thermal stability. In Section 4, we study the radiation re¬ 
leased from disk-like outflows with energy deposition at the 
inner edge and develop a semi-analytic model for luminosi¬ 
ties and effective temperatures. In Section 5, we study the 
dust formation in the outflow. In Section 6, we summarize 
the observable signatures of L 2 mass loss and compare them 
to known red transients. In Section 7, we conclude by dis¬ 
cussing extensions and other possible applications of the de¬ 
veloped framework. 


2 NUMERICAL SETUP 


2.1 Smoothed Particle Hydrodynamics 


We study the structure of mass loss from the L 2 point of 
binary stars using smoothed particle hydrodynamics (SPH). 
The advantages of SPH for this problem are the ability to 
handle anisotropic geometry, large density contrasts, and an 
accurate treatment of nearly-ballistic motion of the particles. 
We employ the conservative formulation of SPH with spa¬ 
tially varying smoothing lengths (Price & Monaghan 2007) 
and perform the calculations in a Cartesian inertial refer¬ 
ence frame with the barycenter of the binary star as the ori¬ 
gin. Our SPH particles move in a prescribed external time- 
varying gravitational field of the binary star in the center 
and we assume that the mutual gravitational interactions 
between the particles are negligible. The evolution of the 
velocity Vj of particle j is 


^ = y, 

d t -4-J 


y.Q.j pfcij 


VW ji + 


rr/; // ; I (//V Wj ; + ^binary. 


( 1 ) 


where m,. p,-, and P , are the mass, density, and pressure of 
particle i, respectively, D; takes into account time-varying 
smoothing lengths, Wji = (W/ + W,)/ 2, where W, = W(rji,hj) is 
the standard cubic smoothing kernel for smoothing length li{ 
and the particle distance is rji = |r,y| = |ry-r,j, where r = ( x,y,z ) 
is the position vector. Following Price & Monaghan (2007), 
we simultaneously solve for hj , pj, and Qj by demanding 
that the mass within the smoothing volume stays constant. 
In each time step, we first predict an initial guess of hj, 
and then refine the value using Newton-Raphson iteration. 
For rare cases when this procedure does not converge, which 
typically occurs for particles immediately after injection at 
L 2 , we use bisection. 

The second sum in Equation (1) represents the artificial 
viscosity to resolve shocks, in particular we use the standard 
SPH viscosity of Monaghan & Gingold (1983), 

( -asPHCsjifti+PsPH/4 n 

%= -- v Jrrji<0, ( 2 ) 

( 0 otherwise. 


where Vj, = Vj - v, is the velocity difference of particles j and 
i, csji = (csj + csj )/2 is the mean adiabatic sound speed, pji = 
(Pj+Pi)/ 2, and 


Pji '■ 


hjiVji- 


r : . + (5sph^ 2 ' 

j, SCH J, 


( 3 ) 


where hjj = (hj + hf )/2 is the mean smoothing length, and 
ffSPH = L0, /?sph = 2.0, and 5sph = 0.01 are coefficients. We 
find that reducing asPH and /?sph by a factor of 10 leads to 
undesired particle interpenetration for our runs. The coeffi¬ 
cient fjj = (fj + f )/2 is the Balsara (1995) “switch”, 


_jV-vj_ 

|V • v| + |V x v| + I0~ 4 csj/hj ’ 


( 4 ) 


which limits the artificial viscosity in shear flows. We find 
that setting fj, = 1 does not noticeably change the results in 
this work, because there is no substantial shear present in 
our setup. We calculate the velocity derivatives for particle 
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j in Equation (4) as 


Z nij _ , , 

jVi-VWji, (5) 

Vxv = J^—ViXVWji. (6) 


We experimented with other possible expressions, but did 
not find any noteworthy differences. 

The last term in Equation (1) is the gravitational accel¬ 
eration due to the central binary 

^binary = ~GM\ *~ - GM 2 ' ~^ , (7) 


where M\, Mi, . and r „2 are the masses and position 
vectors of the components of the binary with total mass M 
and mass ratio q = M 1 /M 2 , M\ < Mi - The components of the 
binary move on a circular Keplerian orbit counter-clockwise 
in the x — y plane with semi-major axis a = |r*q — I and 
orbital period P. The SPH particles do not affect the orbital 
properties of the binary. 

The internal energy Uj of particle j is evolved according 
to 


d Uj 

—— = "gas,2 + "vise ,j "dill,/ + "cool,/ + "heat,/, 


( 8 ) 


where 


gas,7 ' 


pjClj 


E' 


Vji-VWj. 


( 9 ) 


is the standard SPH term describing adiabatic expansion 
and contraction. The viscous heating " v iscJ is calculated as 


'• ^' "lifji 11 v ji ■ V Wj, 


2 x—i 


( 10 ) 


The remaining terms in Equation (8) capture the radia¬ 
tive processes in the outflow. Although radiation transport 
has been previously included in SPH using implicit meth¬ 
ods for flux-limited diffusion (e.g. Whitehouse & Bate 2004; 
Whitehouse et al. 2005), our goal in this work is to estimate, 
for the first time, the dependence of the radiative properties 
of the L 2 mass loss on the parameters of the binary star 
and the mass loss rate. To these ends, we make several ap¬ 
proximations in our treatment of radiative effects focused on 
making the radiative terms as local as possible and hence 
easily evaluated explicitly within the SPH formalism. As¬ 
suming that the radiation and the gas are in equilibrium, 
we split the radiative transport into three independent pro¬ 
cesses: flux-limited radiative diffusion "diff, which only re¬ 
distributes the energy between the particles predominantly 
in the radial direction, radiative cooling w coo i in the vertical 
direction, and the irradiation by the central binary "heat- Ac¬ 
curate predictions of the emergent radiation should be done 
in the future using more sophisticated methods. 

The flux-limited radiative diffusion term "diff.j 
(Bodenheimer et al. 1990; Cleary & Monaghan 1999; 
Mayer et al. 2007; Forgan et al. 2009) is evaluated as 


"diff.y — 




4 /??/ kjki 
P jPi kj kj 


(Tj-Ti)- 


V W,-i 


( 11 ) 


where 7) is the temperature of particle i, and k ,■ is the thermal 
conductivity coefficient described in Appendix Al. Equa¬ 
tion (11) models the diffusion of the radiation as a con¬ 
duction process by interpolating between the optically-thick 


limit producing radiative diffusion and the optically-thin 
limit, when "diff,; vanishes. Since the energy exchange be¬ 
tween particles is pair-wise symmetric, equation (11) con¬ 
serves energy and describes the redistribution of energy 
within the disk-like outflow. The energy losses due to ra¬ 
diative cooling, which occurs primarily perpendicular to the 
orbital plane, must be included separately. 

Stamatellos et al. (2007) and Forgan et al. (2009) ap¬ 
proximate radiative cooling and heating rates for SPH par¬ 
ticles using local values of density, temperature, opacity, and 
gravitational potential. Their work was adapted to more 
anisotropic configurations by replacing the gravitational po¬ 
tential with the gas pressure scale height by Lombardi et al. 
(2015). Neither of these prescriptions are immediately ap¬ 
plicable to L 2 outflows, because our SPH particles are not 
self-gravitating and their final trajectories do not deviate 
dramatically from the ballistic approximation (the outflow 
is not supported by the thermal pressure). Nonetheless, 
inspired by Mayer et al. (2007), Stamatellos et al. (2007), 
Forgan et al. (2009), and Lombardi et al. (2015), we assume 
that the radiative cooling occurs primarily in the direction 
perpendicular to the orbital plane and estimate the energy 
loss rate per particle as 


o- SB T? 


"cool J — ' 


^ZjTzJ + Kj 1 


( 12 ) 


where and t z j are surface density and optical depth 
from the position of particle j out of and perpendicular to 
the orbital plane. This prescription interpolates between the 
optically-thin cooling rate of kjctsbT^ when t z < k 1, to a sup¬ 
pressed cooling rate as t z » 1. In the optically-thick limit, the 
energy is removed approximately on the diffusion timescale 
out of the outflow, " coo i « "/fdiff (Stamatellos et al. 2007). 
Although this suppression of radiative cooling is not very so¬ 
phisticated, it has been successfully used in SPH simulations 
(Stamatellos et al. 2007; Forgan et al. 2009; Lombardi et al. 
2015). However, a more realistic treatment of the radiation is 
definitely warranted. Note that the outflow can be optically- 
thick in the radial direction, yet optically-thin in the vertical 
direction, where most of the cooling occurs. We set u coo \j = 0 
for T < 200 K to prevent particles cooling off the equation of 
state grid. 

In a similar manner to " coo |j, we include heating due to 
irradiation by the central binary as 


"heat, j — 


L, 1 

4nrj-L r jT r j+Kj v 


(13) 


where j and T r j are surface density and optical depth in 
the radial direction. The central binary is treated as a point 
source at the origin radiating isotropically with a bolometric 
luminosity of L* = 4na 2 crsBTt■ The form of Equations (12- 
13) guarantees that in the optically-thin limit the tempera¬ 
ture of a stationary particle converges to the local equilib¬ 
rium temperature given by 



The opacities Kj in Equations (12-13) are technically Planck- 
ian means (Stamatellos et al. 2007; Forgan et al. 2009), but 
for the sake of simplicity we instead use the Rosseland means 
throughout this work. Future refinements in this area should 
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be accompanied by improvements in the radiation trans¬ 
port algorithm. More details of how we calculate the ra¬ 
diative cooling and irradiation heating are provided in Ap¬ 
pendix A2. 

Because of the compact support of the smoothing ker¬ 
nel W, the summations in Equations (1-11) are performed 
only for a modest number of nearby particles. To find the 
nearest neighbors for our particles, we employ ^-dimensional 
trees (Press et al. 2007). The SPH equations are integrated 
using the leapfrog method (Springel et al. 2001). The global 
timestep At is chosen to be the minimum over all particles of 
three timescales: (1) the Courant-Friedrichs-Lewy condition 
(Bate et al. 1995), 


A f CFL = 


_ 03h _ 

Cs + /l|V-v| + 1.2(»sPHdS +/?SPHf ! |V-v|)’ 


(15) 


where the last term of the denominator is included only for 
V ■ v < 0, (2) the force condition 


A f F = 0.3 


h 


|dv/dr| 


(16) 


and (3) the condition that the internal energy should change 
by less than 10% in one timestep, 


At u =0.1 


(17) 


We have checked that our results are not sensitive to modest 
increases in the prefactors in Equation (16-17), allowing us 
to adopt these less computationally expensive values. 

Since the studied problem does not necessitate the cal¬ 
culation of the gravitational forces between the particles, but 
requires frequent injection and removal of particles, query¬ 
ing microphysical tables, and evaluation of global quanti¬ 
ties such as optical depths, we decided to write a new cus¬ 
tom SPH code in the C/C++ language and parallelize using 
OpenMP. In addition to the ^-dimensional tree routines of 
(Press et al. 2007), we use the GNU Scientific Library for 
root finding, sorting, and random number generation. The 
runtime of the application is dominated by the evaluation of 
the microphysics, and the global quantities such as optical 
depths, which is done efficiently in parallel. Construction of 
the particle tree and adding new particles are implemented 
serially, but this requires only very little computing power. 
Our code is explicit, which can lead to a very short, but not 
prohibitive, timestep in situations when the radiative diffu¬ 
sion dominates the energy equation. Typically, we require 
several tens or hundreds CPU hours to reach steady-state 
for optically-thin configurations, but some of our optically- 
thick runs took up to about 3000 CPU hours on a 20-core 
machine. 


2.2 Luminosity and Emission Temperature 

We calculate the total radiative luminosity of the binary 
ejecta as a sum over all particles in the simulation, 

^ ~ 'y \ m j\ucoo\j\- (is) 

i 

The effective temperature of the emission, 7^, is estimated 
from the average of the effective temperatures of individual 
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Figure 1. Adiabatic index T as a function p and T for our equa¬ 
tion of state. Important state changes are labeled in the figure. 


particles, assuming radiative equilibrium in a grey atmo¬ 
sphere and weighted by their radiative energy loss rate, 


“ L^ mj]UcooX ’K z , j + T 


T 4 

j 


eff ' 


(19) 


We emphasize that our treatment of radiative losses and 
its descriptors L and T e g are only approximate. Obtaining 
more accurate predictions of the luminosity and spectral en¬ 
ergy distributions of mass losing binaries in future work will 
require using more sophisticated methods. 


2.3 Equation of state 

We use the equation of state (EOS) of Tomida et al. (2013, 
2015), which assumes equilibrium abundances of ionized, 
neutral and molecular hydrogen including the rotational 
and vibrational degrees of freedom, and neutral, singly- and 
doubly-ionized helium. The composition is X = 0.7, Y = 0.3, 
and the ratio of ortho- to para-H 2 is assumed to be in 
thermodynamical equilibrium. The latent heats for all pro¬ 
cesses are consistently included. The EOS is given on a 
densely-spaced logarithmic grid spanning -22 < logp < 1.1 
and 0.2 < log T < 6. We modify the gas EOS by adding the 
contribution of radiation to the pressure and energy. How¬ 
ever, we also need to modify the adiabatic sound speed cs , 
which enters the calculation of artificial viscosity and time 
stepping conditions. Based on time-steady Euler equations, 
Pejcha & Thompson (2012, Appendix A) expressed cy as a 
function of only derivatives of P and u with respect to ei¬ 
ther p or T, while keeping the other quantities fixed. These 
thermodynamic derivatives can be easily evaluated on the 
gas EOS grid and we verify that the resulting gas cs is re¬ 
produced with better than 1% accuracy. We then repeat the 
procedure with radiation contribution included in P and u 
to obtain the final cs . Finally, we remap the EOS grid from 
the p-T to p-u plane to reduce the runtime computational 
overhead. At runtime, we obtain the thermodynamic quan¬ 
tities of interest by bi-log-linear interpolation. 

In Figure 1, we summarize the basic features of the fi¬ 
nal EOS using the adiabatic index T = ( dlnP/dlnp)s = pc| /P. 
At low densities, the EOS is dominated by radiation with 
F = 4/3, while at higher densities the gas physics results in 
1.0 < T < 5/3. Specifically, for logT < 1.5, the EOS behaves 


MNRAS 000, 000-000 (2015) 













6 Pejcha et al. 


as a mono-atomic gas with T ~ 5/3, because the temper¬ 
ature is not high enough to excite rotational degrees of Hb 
molecules. For 1.5 < log 7" < 2.0, the conversion between para- 
and ortho-H 2 reduces T. For 2.0 < logT < 3.0, we find T » 7/5 
corresponding to an ideal diatomic gas. For logT > 3.0, the 
latent heats associated with IF dissociation, and H and He 
ionization result in T « 1 for relatively large regions of the 
parameter space. The practical result is that if the energy 
input due to shock heating is comparable to the binding 
energy of H 2 (2.26 eV per baryon) or the ionization energy 
of hydrogen (13.6eV), the temperatures will pile up between 
1500 and 3000 K or between about 8000 and 12000 K, respec¬ 
tively. This is important for trying to analytically estimate 
the shock temperatures, as we discuss in Section 4. 


2.4 Opacities 

Opacities enter our calculation through the diffusion co¬ 
efficients (Eqs. [11], [Al]), the estimates of radiative 
cooling (Eqs. [12], [A6]), and irradiation by the cen¬ 
tral star (Eq. [13]). We employ the solar-metallicity Z = 
0.0189 opacity tables compiled in the stellar evolution 
code MESA version 7503 (Paxton et al. 2011, 2013) by 
modifying the tutorial example sample_kap. f in the kap 
module. Using this code we produce two opacity tables 
for the two possible choices of low-temperature opacity 
source by changing kappa_low_T_pref ix option. The op¬ 
tion lowT_f a05_gn93 produces low-temperature opacities of 
Ferguson et al. (2005), which include the effects of molecules 
and grains. The option lowT_Freedmanll yields grain-free 
molecular opacities of Freedman et al. (2008) with up¬ 
dates to the molecular hydrogen opacity (Frommhold et al. 
2010) and ammonia opacity (Yurchenko et al. 2011). Both 
low-temperature opacity tables are blended with OPAL 
(Iglesias & Rogers 1993, 1996) at high temperatures, as de¬ 
scribed in the MESA instrument papers (Paxton et al. 2011, 
2013). Opacities outside of the boundaries of the input ta¬ 
bles are set to the values at the boundary. Although this 
potentially introduces bias in our calculations, it is prob¬ 
ably rather small, because most of the radiation emission 
occurs at temperatures and densities covered by the tables. 

The final opacity tables are visualized in Figure 2. 
For both tables, the opacity steeply rises for logT > 3.7 
due to H~ opacity. The opacities are very low for 3.1 < 
logT < 3.6 and 3.3 < logT £ 3.8 for lowT_f a05_gn93 and 
lowT_Freedmanll, respectively. However, the tables signifi¬ 
cantly differ for logT <3.1, where the presence of dust grains 
in the lowT_f a05_gn93 table results in relatively high opac¬ 
ities. Conversely, the grain-free table lowT_Freedmanll ex¬ 
hibits very low opacities for low temperatures. Since the dust 
formation does not depend exclusively on the temperature, 
but also on other properties of the outflow (Section 5), we 
employ both opacity tables and compare the results. 


2.5 Initial and boundary conditions 

Particles are injected around the L 2 point with initial y and 
z coordinates drawn randomly from a Gaussian distribution 
of width sa, where £<K 1 estimates the typical size of the L 2 


region and is defined following Shu et al. (1979) as 



where cj = JkT/m p is the isothermal sound speed, m p is 
the proton mass, and v c j rc = x/GM/a is the binary orbital 
velocity. The corresponding x coordinate for each particle 
is fixed by requiring that it lies on a sphere centered on 
M\ and intersecting the L 2 point. We adopt a fiducial value 
of e = 0.05 (Eq. [20]) and investigate the sensitivity of our 
results in Appendix B. The radial coordinate of the L 2 point, 
tl 2 = r^a, is determined by numerically solving the usual 
equation (e.g. Shu et al. 1979). The initial velocities of the 
particles are set to the corotation velocity at the L 2 point, 
vl 2 = 2nr Ll /P. 

A particle is only injected into the full simulation if its 
purely ballistic motion keeps it exterior to L 2 for a moderate 
initial ‘trial’ period. We experimented with different meth¬ 
ods of injection, such as placing particles directly on the 
Roche surface, starting the stream slightly offset from L 2 
or imparting particles with an initial radial velocity kick of 
the order of svl, • The results in each case were nearly iden¬ 
tical, indicating that the specific injection procedure is not 
important. For a given particle injection rate N and target 
mass loss rate M, we determine the total mass of all active 
particles at each time step. If the actual number of particles 
is lower than the target, additional particles are added with 
constant mass M/N. Although the number of added particles 
is slightly different in each time step depending on the num¬ 
ber of particles removed from the simulation, this procedure 
guarantees that we achieve the desired global time-averaged 
M. Typically, we choose N = 1000/T as our low-resolution 
calculation and explore the sensitivity of our results to N in 
Appendix B. 

Even if the ballistic trajectory of a particle would other¬ 
wise take it outside of L 2 , pressure and viscous forces from 
neighboring particles might still drive it back inside. Such 
particles would frenetically orbit the two point masses inside 
L 2 , significantly reducing the global integration timestep and 
occasionally resulting in the violent ejection of the particle. 
We account for this clearly unphysical behavior by artifi¬ 
cially removing any particles that reside closer to the co¬ 
ordinate origin or to M\ than their injection position (usu¬ 
ally such removals happen within a few timesteps after in¬ 
jection). Typically, roughly one third of the initial parti¬ 
cles survive to participate in the global dynamics. The suc¬ 
cess fraction of particle injection increases if the injection 
point is moved slightly outside of L 2 , although we find that 
such modest changes do not significantly alter the global 
outcome of the calculations. Our choice of an inflow in¬ 
ner boundary condition motivates our focus on binaries 
with 0.064 < q < 0.78, which produce unbound outflows (Sec¬ 
tion 3.1, (Section 3.1 and Shu et al. 1979). 

After particles are injected, we fix their temperatures 
artificially at the equilibrium temperature set by binary ir¬ 
radiation (eq. 14) until the radial coordinate of the parti¬ 
cle exceeds 1.4rL 2 , corresponding to a small fraction of the 
first spiral wrapping. This approach removes transient ef¬ 
fects that occur because some particles are injected super¬ 
sonically with respect to others, in particular by mitigating 
shock heating within the stream very close to L 2 . Although 
such internal shocks are not obviously unphysical, our global 
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2.0 2.5 3.0 3.5 4.0 4.5 5.0 2.0 2.5 3.0 3.5 4.0 4.5 5.0 

log T log T 


Figure 2. Opacity tables as a function of temperature and density. Left panel shows low-temperature opacities from Ferguson et al. 
(2005) (option lowT_f a05_gn93 in MESA), while the right panel is using opacities from Freedman et al. (2008) (option lowT_Freedmanll). 
Outside of the boundaries of the input tables the opacities are set to the value at the boundary, which leads to stripes at low temperatures. 


results are not sensitive to the precise inner temperature evo¬ 
lution. We have experimented with keeping the particles at 
a constant temperature (instead of T eq ), or not masking this 
transient at all, finding no lasting influence on the structure 
of the stream beyond the first spiral wrapping. 

3 STRUCTURE OF THE L 2 OUTFLOW 
3.1 Mass loss from Li 

Here, we summarize results on the semi-analytic treatment 
of mass loss from the outer Lagrange point L 2 , which is based 
on the work of Shu et al. (1979). Mass loss begins once the 
common surface of the synchronously rotating contact bi¬ 
nary star comes into contact with L 2 (Kuiper 1941). Well in¬ 
side of this surface, gas is supported in hydrostationary equi¬ 
librium by a balance between pressure and effective gravity 
(combination of gravitational and centrifugal accelerations). 
At L 2 the confining force vanishes, allowing gas to transition 
to a hypersonic outflow. As argued by Lubow & Shu (1975) 
for the Li point and Shu et al. (1979) for the L 2 point, the 
gas-emitting region around L 2 has a characterstic size ea, 
with the transition to supersonic motion occuring on a sim¬ 
ilar spatial scale. The value of e depends only weakly on 
the parameters of the binary (Eq. [20]). Its small value in¬ 
dicates that, away from the immediate neighborhood of L 2 , 
the thermal pressure of the star does not significantly affect 
the kinematics of the outflow (ballistic motion). 

As found by Shu et al. (1979), only binaries with 0.064 < 
q < 0.78 torque the gas enough to guarantee ejection to in¬ 
finity. We confirm this result by following the motion of a 
number of ballistic particles randomly injected near L 2 from 
a neighborhood of size ea. The particles are assigned random 
velocity perturbations of the order of svl 2 , where vl, is the 
corotation velocity at L 2 . Figure 3 shows the fraction of par¬ 
ticles that successfully escape to infinity. We find that > 50 
per cent of the particles are unbound for q within the limits 
defined by Shu et al. (1979). Figure 3 also shows the mean 
asymptotic velocity Voo of the unbound particles, which are 
typically between 0.15-0.25 of the binary escape velocity, 
v eS c = V2 GM/a. Our results for Voo also closely match those 



Figure 3. Fraction of unbound particles ballistically ejected from 
the £ vicinity of L 2 (black solid line and circles) and their mean 
asymptotic velocity (blue stars) in the units of - \J2GMIa as a func¬ 
tion of the mass ratio of the binary. Dashed blue line shows asymp¬ 
totic velocities from Table 2 of Shu et al. (1979). The vertical dot¬ 
ted lines indicates the region, where Shu et al. (1979) found the 
stream escapes, 0.064 <q< 0.78. 


of Shu et al. (1979), except for values of q with only a small 
fraction of unbound particles. Since the ballistic problem is 
scale-free, these results are valid for all binaries. The rela¬ 
tive insensitivity of Voo to q further suggests that measuring 
the L 2 outflow velocity (from spectroscopy of a brightening 
merger event, for example) will constrain v es c and thus the 
total mass and semi-major axis of the binary. 

In the rest of the paper we focus on binary mass ratios 
0.064 < q < 0.78 which produce unbound L 2 outflows. For 
binaries with q < 0.064 or q > 0.78, we confirm the finding of 
Shu et al. (1979) that gas piles up at a finite outer radius 
and forms a bound mass-loss ring instead of an unbound 
outflow. 
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Figure 4. Structure of the outflow from Lo for a binary with M = 1.725M Q , q = 0.15, a = 0.03 AU, and M = 10 3 /V/ Q yr 1 . The top left panel 
shows the surface density, while the remaining panels show the temperature for progressively more added physics. The potential surface 
containing L 2 is indicated with a red line near the center. 

3.2 Density and temperature structure 

In Figure 4, we show the typical density and temperature 
structure of the outflow from our radiation-hydrodynamics 
calculations. Initially, the stream of gas leaving L 2 fol¬ 
lows the spiral ballistic trajectory (Shu et al. 1979), while 
it expands and cools. The difference of gravitational forces 
across the stream contributes to the stream expansion. If 
the stream is optically-thin, then irradiation by the cen¬ 
tral binary compensates for the radiative cooling and adia¬ 
batic losses and the temperature converges to T eq (r) (eq. 14). 

For optically-thick streams, both irradiation and cooling are 
comparatively unimportant and adiabatic losses dominate. 

Although the stream expands significantly, the temperature 
drop is typically smaller than expected from a simple poly¬ 
tropic EOS due to the release of latent heat from H 2 forma¬ 
tion. 

When the spiral wraps around the binary » 1.5 times 
and reaches the radial distance of r c ~ 8 a, the outer edge of 


the inner spiral encounters the inner edge of the previous 
wrapping of the spiral. 1 The associated velocity difference 
produces a shock, which heats part of the outflow. The ther¬ 
modynamic properties at r > r c are set by the vertical opti¬ 
cal depth at r c and the relative importance of the individual 
terms in the energy equation. We illustrate this in Figure 4 
by showing how the outflow temperature and density struc¬ 
ture changes as the individual terms in Equation (8) are 
switched on one by one. 

If the radiative diffusion and cooling are unimportant, 
the shock-heated material is advected outward forming a spi¬ 
ral (top right panel of Figure 4). Typical peak temperatures 
are ~ 10 4 K, depending on binary properties. If instead 


1 Note that r c depends on the initial spread of the stream and 
its temperature: if the stream was initially cold, e = 0, no colli¬ 
sion would occur between spiral wraps. We can thus approximate 
r c /a ~ . 
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Figure 5. Outflow density as a function of radius r for the same 
binary as in Figure 4, but for a range of mass loss rates M. Black 
lines show the median of the particle distribution, while colored 
bands show the 5-95% percentile range. For r> 10a the density 
profile obeys p oc r~ 2 . 


dominates the energy evolution and radiation cooling is neg¬ 
ligible, then the shock-dissipated energy is efficiently redis¬ 
tributed within the outflow and the spiral structure is erased 
(bottom left panel of Figure 4). If the outflow is optically- 
thin in the vertical direction, then radiative cooling domi¬ 
nates, and the shocked gas cools almost immediately at the 
position where the two spiral wrappings first intersect (bot¬ 
tom right panel of Figure 4). If the outflow is optically-thick, 
then the heated gas needs to be advected to larger radii be¬ 
fore it can efficiently radiate and part of the thermal energy 
is then lost to adiabatic expansion. 

In Figure 5, we show the density profiles as a function 
of the radial distance r from the binary center. Initially, 
at small r the density profile decreases steeply, p oc r~ 4 . For 
r > r c ~ 8a, the density profile then flattens and behaves ap¬ 
proximately as pocr 2 . This asymptotic behavior is useful 
in constructing a semi-analytic model in Section 4. There 
is no noticeable increase of density around r c . There are no 
remarkable differences in the density profiles for different bi¬ 
nary configurations and mass-loss rates, except for its overall 
normalization. 

In Figure 6, we show the vertical distribution of particles 
as a function of distance in the orbital plane g = -\jx 2 + y 2 and 
compare this to the analytic prediction of Shu et al. (1979, 
Eq. 61) 2 , 


Z°c 


( g/a)e yjlnlnig/aft. 


( 21 ) 


Agreement between our numerical results and Equation (21) 
is good if the constant prefactor is approximately unity. Note 
that the thickness of the outflow increases linearly with the 
relative size of the L 2 region s (Eq. [20]) and faster than 
linearly with g. The outflow thickness is not found to depend 
sensitively on M or q. For the fiducial e = 0.05, the typical 


2 We rewrite the result of Shu et al. (1979) by assuming that 
the particles move radially with a constant velocity. This is good 
approximation for r>r c ~ 8a, as we in Section 3.3. 


Figure 6. Vertical distribution of particles as a function of the 
distance g = fx 2 + v- in the orbital plane, calculated for the same 
binary as in Figure 4. The shaded area encompasses 5% and 95% 
of the particles at a given distance in the orbital plane. Dashed 
lines show the expected growth of the thickness for free isothermal 
expansion in the vertical direction given by Eq. (21) (Shu et al. 
1979). 


half-opening angles at r c are between 5 and 10 degrees and 
the typical vertical thickness is ~ a. 

Finally, we investigated the effects of irradiation by the 
central binary and the alternative opacity tables on the 
results. Although irradiation heating primarily affects the 
temperature structure of optically-thin spiral streams, these 
changes are typically confined to the innermost spiral wrap¬ 
ping and do not result in qualitative changes in the spiral 
merging and associated energy release. Optically-thick out¬ 
flows are found to reprocess up to about 10 per cent of the 
binary luminosity. We also compared our results for a series 
of calculations employing low-temperature opacities from 
lowT_fa05_gn93 to those using lowT_Freedmanll. We find 
differences which are generally less than an order of magni¬ 
tude, because most of the emission occurs under conditions 
of density and temperature where both tables provide sim¬ 
ilar results. The grain-free lowT_Freedmanll opacities gen¬ 
erally produce higher luminosities and slightly higher effec¬ 
tive temperatures, because the lack of grains allows trapped 
radiation to escape earlier. Since the conditions for dust for¬ 
mation are generally favorable in the outflow (Sec. 5), we 
hereafter primarily employ the lowT_f a05_gn93 opacity ta¬ 
ble, which includes condensates. 


3.3 Source of shock energy 

We now address the origin of the energy dissipated by the 
shock. Since the gas is moving highly supersonically, the 
kinetic energy dominates the total energy budget. In Fig¬ 
ure 7, we show the velocity distribution of particles in one 
simulation snapshot and compare them to a range of purely 
ballistic trajectories with similar initial conditions. From the 
total velocity in the left panel, we see that the distribution of 
SPH particles, specifically the spread in velocity at a given 
radius, is initially nearly identical to ballistic particles. The 
particles are torqued and accelerated by the time-changing 
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Figure 7. Velocity distributions of particles as a function of radius, calculated for the same binary parameters as in Figure 4. We 
show the total velocity (left panel), radial component |v||| = v r/r (middle panel) and the tangential component [v L | = [v-vy (right panel). 
Black lines show the median of the particle distribution, while green shading marks the 5 - 95% percentile range of the distribution. 
For comparison we show purely ballistic particle trajectories with similar initial conditions (pink hatched regions), and the local escape 
velocity V2 CfMjr (red dashed line). 


gravitational field of the binary, as evidenced by the growing 
disparity between the particle velocities and the local escape 
velocity V2GM/?\ After particles are shocked at r c ~ 8 a, the 
velocity distribution becomes narrower, implying that faster 
moving particles transfer their energy to the slower moving 
particles. The velocity change Av is about 10 to 20 per cent of 
the escape velocity at r c ~ 8 a. In the middle and right panels 
of Figure 7, we show the radial and tangential components 
of the velocity. At r c , the radial component dominates over 
tangential. There is little change in the tangential velocity 
distribution at r c and the distribution is narrow at all r. 
This implies that shock dissipation energy comes from the 
homogenization of the radial velocities of particles and that 
there is very little shear. This also explains why the simu¬ 
lation results do not depend on whether the Balsara (1995) 
artificial viscosity “switch” is turned on or off. 


3.4 Instabilities in the outflow 

The outflow can be unstable to the gravitational formation 
of clumps, which could potentially violate our assumption 
of negligible gravitational interactions between the parti¬ 
cles. The required condition is that the gravitational collapse 
time fg lav ~ (Gp)~’/ 2 is shorter than the expansion timescale 
f ex p ~ r/voo . Assuming that M ~ r~pv oo (Fig. 5), the stability 
condition is 

M \ 1/2 / vpo 
O.lMoyr^j \100kms~ 1 

The highest luminosities of the binary outflow investigated 
in Section 6 require not only high M but also high v e sc, 
and Equation (22) implies that even these outflows are sta¬ 
ble. This also justifies our assumption that the self-gravity 
among the particles is not important. The stability condi¬ 
tion roughly holds at r c even if we replace the expansion 
timescale with the sound crossing timescale. 

Alternatively, the outflow might form clumps by becom¬ 
ing thermally unstable. The condition for thermal instabil¬ 
ity roughly requires that the cooling rate increases along the 


\— 3/2 

« 1 . ( 22 ) 



thermodynamic trajectory of the gas parcel. This can hap¬ 
pen during dust formation, when a small decrease of tem¬ 
perature leads to the increase of k, which can yield higher 
cooling rate, more dust formation, and thus further increase 
in k. In the context of multi-phase interstellar or intergalac- 
tic medium, conduction suppresses the formation of clumps 
smaller than the Field length 

(23) 

y P u cool 

where k is the conduction coefficient (e.g. Field 1965; 
Begelman & McKee 1990). In this work, the energy in the 
outflow is redistributed by radiative diffusion, which is mod¬ 
eled as conduction in our SPH code (App. Al). The binary 
outflow becomes thermally unstable if /If <k r and the clumps 
fit within a causally connected region. We evaluated dp in 
several runs covering various values of M and v esc and find 
that unless the radiative diffusion dominates the evolution, 
parts of the outflow are thermally unstable with Ap/r be¬ 
tween 0.01 and 0.1. Our simulations do not have sufficient 
resolution to demonstrate whether the clumping actually oc¬ 
curs or to study the properties of the clumps, but this topic 
deserves further attention. Clumping of the outflow might 
be responsible for the small quasi-periodic variations in the 
rising light curve of V1309 Sco (Pejcha 2014). 


3.5 Backreaction on the binary orbit 


Gas leaving the binary carries away the angular momentum 
of the L 2 point, j = 2nM{rp 2 a) 2 1P. If the orbit remains cir¬ 
cular, the orbital period of the binary shrinks as 



(24) 


where 3K strongly depends on the mass ratio of the binary 
and on how does q change throughout the mass loss (e.g. 
Pribulla 1998). 

The gas leaves to infinity relatively “cleanly” in the sense 
that it interacts with the central binary only through the 
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time-changing gravitational field and not directly through 
collisions and shocks around the stars. To assess the magni¬ 
tude any potential additional angular momentum loss due 
to gravitational torquing, we calculate the total angular mo¬ 
mentum of the outflow, J = Tjj m j r j x Vj, and compare its 
value to the sum of angular momenta of particles at the 
time of their injection J j n j t . We find that \J\ very soon asymp¬ 
totes to values about 10 to 18% higher than |/j n i t |, depending 
on q. This is a relatively small correction to 3K, especially 
considering the much greater influence of uncertainties in q 
on JK (Pribulla 1998). For binaries forming a mass-loss ring 
(q < 0.064 or q > 0.78), the angular momentum extraction 
might be much more efficient. 


4 RADIATIVE PROPERTIES AND 

SEMI-ANALYTIC MODEL 

4.1 Results from simulations 

In Figure 8, we show the asymptotic luminosities and ef¬ 
fective temperatures obtained by simulating binaries with 
constant M long enough for the outflow to establish a steady- 
state. The asymptotic luminosity is approximately linearly 
proportional to M , while binaries with higher escape veloc¬ 
ities are also more luminous. For the range of parameters 
considered in Figure 8, luminosities reach up to 10 6 L o - Note 
that the steady-state assumption may not be fully appropri¬ 
ate to describe some of the red transients, a point we discuss 
further in Section 6. 

The right panel of Figure 8 shows the effective tem¬ 
perature of the emission. For the range of binary properties 
shown, 500 £ T e fj < 6000 K. Although depends only very 
weakly on the mass-loss rate, a clear trend emerges in that 
binaries with higher escape velocities v es c produce higher 
temperature emission. This agrees with the physical picture 
presented in Section 3.3 that the shock dissipates a fixed 
fraction of the outflow kinetic energy oc v 2 sc (Section 3.1). 
The luminosities and effective temperatures derived from 
our simulations are similar to those of the observed red tran¬ 
sients (Sec. 6). 

In Figure 9, we show the asymptotic luminosity and 
effective temperature as a function of mass ratio, when all 
other parameters are held fixed. We see that both L and 
r eff exhibit minimum at q ~ 0.3, and increase on either side. 
This behavior can be understood by analyzing the motions 
of ballistic particles (Fig. 3), which show that the fraction of 
particles escaping to infinity decreases when q is either large 
or small. As the efficiency of particles escaping to infinity 
decreases, the collisions become more efficient, resulting in a 
higher luminosity. Since the luminosity is released at roughly 
the same distance from the binary, T e ff follows the same 
trend as L. In all cases, the luminosity varies by less than a 
factor of 2 and hence the dependence on q can be neglected 
in analytic estimates. 

In Figure 10, we show the dependence of the luminosity 
and effective temperature on the total mass of a binary with 
a and q held fixed. As expected, L and T e ff increase with 
M, because the binary has higher v esc and therefore the ve¬ 
locities of the collisions are higher. Equations (25) and (29) 
suggest that L oc M, roughly similarly to what is observed in 
Figure 10. 


A striking feature is the sharp drop in L at M » 3 M Q . 
We have carefully verified that this feature is not a numer¬ 
ical artifact 3 , but instead represents a transition between 
outflows which efficiently form dust (M < 3 M e ), and those 
in which radiative diffusion dominates (M > 3 M 0 ), keeping 
the outflow isothermal and preventing dust formation. Note 
that the boundary of this transition depends sensitively on 
the semi-major axis, and other binary properties. We inves¬ 
tigate the thermodynamic structure of these outflows and 
the role of dust formation in Section 5. 


4.2 Semi-Analytic model 


The complex range of physics (e.g., EOS, opacities) at 
work and their interplay with radiative diffusion and cool¬ 
ing makes it impossible to construct a fully analytic model 
that would reproduce our simulation results to better than a 
few orders of magnitude. Nevertheless, to highlight the per¬ 
tinent physics and to provide order-of-magnitude estimates 
relevant for comparing with observations, we construct a 
semi-analytic model, which predicts the time evolution of L 
and r eff for a given binary by making use of only the EOS 
and opacity tables. 

For two streams to collide at a radius r c = qa ~ 2qR t 
the velocity dispersion between them must be Av/voo ~ 2?/ -1 , 
where Voo ~ 0.25v es c is the typical asymptotic velocity (Fig. 3) 
giving Voo ~ 50kms _1 (m./rs,) 1 ^ 2 , R„ ~ a/2 is the radius of the 
primary, and // — 8 (Fig. 4), as determined from the sim¬ 
ulations 4 . Here the primary mass and radius are normal¬ 
ized to Solar values, r$ t =R„/5R Q and m* = M/M Q . The gas 
gains energy Ah ~ Av 2 /2 during the collision, which raises 
its temperature. Neglecting H 2 destruction or H ionization, 
the nominal shock temperature for an ideal gas would be 
(Laniers & Cassinelli 1999, p. 360) 


3 o 

T s h “ 77 t Ay2 ~ 3500K 
16 k B 




However, for many binary configurations the latent heat as¬ 
sociated with dissociation or ionization of hydrogen is cru¬ 
cial for determining T s h. In the left panel of Figure 11, we 
show T s h for a range of Av using our EOS and for sev¬ 
eral typical densities. T s h is nearly constant at ~ 2000 to 
45001< for 5 < Av < 20kms _1 due to the latent heats of H 2 
formation and destruction. Similarly, the behavior of T s h 
for 20 < Av < 50kms _1 is affected by hydrogen ionization for 
p > 10~ 10 gcm~ 3 . For lower densities, radiation dominates the 
EOS, resulting in a weaker dependence of T s ^ on Av, as illus¬ 
trated by the p= 10~ 12 gcm -3 case shown in Figure 11. As 
a result of these complications, the correspondence between 
T sh and Av matches Equation (25) only roughly (dotted line 


1 The precise value of M where the drop occurs is found to depend 
slightly on the simulation resolution and the exact procedure used 
to calculate spatial derivatives in our SPH equations. 

4 The parameter q is related to e, because we would expect q —» 00 
for e —* 0, when no collisions between the streams occur. However, 
s describes the thermal content of the gas at the injection point, 
which is then subsequently modified by irradiation, adiabatic and 
radiative cooling, and especially the torqueing by the gravita¬ 
tional field of the binary. As a result, the relation between £ and 
q is only indirect and we keep them separate in the semi-analytic 
model. 
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Figure 8. Asymptotic radiative cooling luminosity (left panel) and the effective temperature (right panel) as a function of (temporally 
fixed) mass loss rate M for simulations which have achieved a steady-state outflow. Different binary parameters (M,a, v eS c) are denoted 
with symbols as described in the legend. Irradiation by the central star was not included so as not to dominate the shock-powered 
emission at low M. In some cases the time needed to reach a steady state exceeds the time over which the binary can lose mass at the 
rate M; we take this constraint into account when discussing applications of our results to red transients in Section 6 . 







Figure 9. Luminosity L (black circles) and effective temperature 
T e ff (red stars) as a function of mass ratio q for the binary from 
Figure 4. The outflow is radiatively cooling in the optically-thin 
regime in all cases. 


in Fig. 11). In the semi-analytical model, we calculate T s h 
directly from Av using our EOS. 

The vertical surface density of the outflow at radius r is 
approximately given by 




M 

InrVoo' 


(26) 


where v TC ~ 0.25v es c (Figs. 3 and 7). The vertical optical depth 
t z ~ "L z k through the outflow at the collision radius r c is 


Mk 
2 nr c v a 


: 3 f ■ n _ 2 K 2 -i Ws.r 1/2 ift- 3 , (27) 

\ 10 z cm z g 1 / 

where M = M_ 3 10 — 3 M 0 yr _1 , and the opacity is evaluated at 
(T s h,p c ), where p c is the density at r c , 



(28) 


Figure 10. Dependence of the luminosity (black circles) and ef¬ 
fective temperature (red stars) on the total binary mass M for 
fixed a = 0.01 AU, <j = 0.15, and M = 10 ~ 3 M 0 yr~*. The jump in curve 
is caused by the transition between diffusion-dominated and dust- 
dominated outflows, as explained later. 


Here we have assumed that the vertical thickness of the out¬ 
flow at r c approximately equals a for e~ 0.05 (Sec. 3.2, Fig. 6, 
and Eq. [21]). Although the outflow spreads in the vertical 
direction, using p c from Equation (28) work well in the semi- 
analytic model. Note that in our actual simulations we cal¬ 
culate 2- and t z more accurately than in Equations (26-27), 
as described in Appendix A2. 


4-2.1 Optically-thin outflow 

If the outflow is optically-thin in the vertical direction, r ZtC < 
1, then its luminosity directly matches that dissipated by the 
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Figure 11. Left : Temperature of the post-shock gas as a function 
of the shock velocity Av « V2Aw, shown for several (fixed) densi¬ 
ties. The ideal gas approximation from Equation (25) are shown 
as dotted lines. Right : Opacity corresponding to the temperatures 
and densities from the left panel. 



Figure 12. Time evolution of the luminosity (solid lines) for sim¬ 
ulations of a binary from Figure 4, shown for different values of the 
(temporally fixed) mass loss rate M as marked. Shown for compar¬ 
ison is the best-fit model using Equation (31) (dashed lines), and 
the predictions of the semi-analytic model (dotted lines, Eqs. [31- 
33]). 


shock, 


Ahin ~ ImAv 2 ~ 5 x 10 34 ergs 1 M _ 3 j ~ j m f rf. (29) 

When H 2 dissociation or H ionization can be neglected, the 
effective temperature of the emission is given by 

T eff ~ ' ~ 2500K M 3/ 3 V /4 r- 3/4 (^) ', (30) 

where cr is the Stefan-Boltzmann constant. As shown in Fig¬ 
ure 4 (bottom right panel), the optically-thin cooling is con¬ 
centrated in a narrow half-ring of radius ~ r c and not the 
full disk assumed in (eq. 30); r e ff nevertheless depends only 
weakly on such geometric factors. 



Figure 13. Ratio of the steady-state luminosity to the optically- 
thin luminosity (Eq. [29]) as a function of optical depth at the 
collision radius r z>c . Different symbols mark series of models with 
constant M, q, and a, but varying M. All calculations were per- 
fomed without the irradiation by the central binary to not bias 
the results due to reprocessing of the central binary radiation. 


4-2.2 Optically-thick outflow 


Predicting the luminosity of an optically-thick outflow with 
T z,c <; 1 is not as straightforward as in the optically-thin case. 
Because an optically-thick outflow cannot immediately ra¬ 
diate energy dissipated by the shock, part of the shock lu¬ 
minosity is lost to PdV work as the outflow expands. The 
details of this process depend on the opacities and EOS, 
which are not simple functions of density and temperature 
(Figs. 1 and 2). We first attemped to predict the asymptotic 
luminosity based on simple assumptions (polytropic EOS, 
constant k), but found several orders-of-magnitude discrep¬ 
ancies between the analytic and numerical results. We there¬ 
fore take a more empirical approach based on a set of our 
simulation results. 

For fixed binary parameters and mass outflow rate M, 
the time evolution of luminosity from our simulations is well 
described by the functional form 


L(t) — /-thick ,max 


1 - exp 


.2 i 


scale ' 


(31) 


where Lthick,max is the asymptotic luminosity achieved at late 
times, and / sca i e is the timescale to approach the steady state. 
For t <K / S cale, Equation (31) gives the expected results that 
the luminosity increases in proportion to the surface area 
of the outflow, Loct 2 . Equation (31) applies also to the 
optically-thin case, assuming we adopt a small value of /scale- 
Figure 12 compares our simulation results for L(t) to equa¬ 
tion (31) for several values of M, spanning the optically-thin 
to optically-thick regimes. We generally End good agreement 
between the simulation results and model fit. 

Figures 13 and 14 show how the parameters in Equa¬ 
tion (31) depend on the key properties of the outflow. 
Figure 13 shows the ratio of the asymptotic luminosity 
L t hick,max to the hypothetical optically-thin luminosity L t hi n 
from Equation (29) as a function of the vertical optical depth 
T z (r c ). For t Zic < 1 we obtain L t hi c k,max/Fthin ~ 1 as expected, 
but for T- c <: 1 the asymptotic luminosity is suppressed with 
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Figure 14. Same as Figure 13 but for the timescale to approach 
the steady state / sca i e . 


increasing r ZjC , approximately as a power-law. This behavior 
is well fit by the functional form, 

f-lhickjiiax ~ f-lliin [ max( I, T z c ) | (32) 

Although the value of the power-law index (-0.2) is not 
obvious from any order-of-magnitude estimates, this is not 
surprising given the complex interplay in our simulations 
between radiative diffusion and cooling, and the properties 
of the EOS and opacities. The scatter about the best-fit 
power-law is smaller than order of magnitude, which indi¬ 
cates that our semi-analytic model captures the pertinent 
processes with good accuracy using Equation (32). Part of 
the scatter comes from the difference between the diffusion- 
dominated and cooling-dominated outflows described in Sec¬ 
tion 5. We emphasize that our fit to Tthick,max( T z,c) factors in 
the size and mass of the binary, as well as the mass-loss rate, 
each of which varies by orders of magnitude across the range 
of simulated stellar systems. 

Figure 14 shows how the characteristic timescale to 
reach steady state f sca i e depends on t Z jC . As with Tthick,maxi 
the f sca i e is approximately constant for t Z jC < 1, while increas¬ 
ing as a power law for r ZiC > 1. This behavior is well fit by 


fscaie ~ — [max(l,r z c)] 0,18 
Voo 


(33) 


across all r ZjC . Although the scatter along Equation (33) is 
again less than an order-of-magnitude, this can introduce 
potentially greater inaccuracies when trying to predict L(t) 
for t < ? sca j e from Equation (31). 

In Figure 12, we show the predicted light curve calcu¬ 
lated from our semi-analytic model ( equations [31] [33]) 
in comparison to the actual simulation results. The agree¬ 
ment is relatively good. For other binary parameters, the 
agreement can be noticeably worse, for example due to time- 
dependent dust formation effects (Sec. 5). Nevertheless, the 
semi-analytic model typically matches the simulations com¬ 
parable to or better than order-of-magnitude for t > t sc ale- 

Finally, we estimate the asymptotic effective tempera¬ 
ture at late times t » f sca i e in the optically-thick case as 


f-thick.i 


[ eff 


47t(v oo fscale)“ cr SB 


(34) 


Figure 15. Comparison of asymptotic effective temperatures 
measured in the simulations and predicted by the semi-analytic 
model (Eq. [34]). The symbols are the same as in Figures 13 and 
14. 



Figure 16. Schematic configuration of mass-losing binary stars, 
where our numerical and semi-analytic results describe the re¬ 
sulting electromagnetic emission. Binary star with angular mo¬ 
mentum vector J and escape velocity v esc = \/2GM/a produces 
an equatorially-concentrated outflow with an asymptotic veloc¬ 
ity Voo x v esc , which is heated at a distance r c from the binary with 
a rate M Av 2 , where Av oc v esc ■ The radiation escapes the outflow 
primarily in the vertical direction. 


In Figure 15, we compare the effective temperatures from the 
simulations to those predicted by Equation (34). The semi- 
analytic model underpredicts the temperatures, especially 
for r eff > 3000 K, an important fact to keep in mind when 
comparing our semi-analytic model to observations (Sec. 6). 


4-2.3 Extensions of the model 

The semi-analytic model presented in this Section can be 
easily extended to a more general class of models relevant 
for the dynamic shrinking of the binary orbit inside of a com¬ 
mon envelope. The semi-analytic model describes a disk-like 
primarily radial outflow with asymptotic velocity Voo, which 
is heated at the inner edge (r ~ r c ) at a rate L^ m (Eq. [29]), as 
shown in Figure 16. The physics of L 2 mass loss determines 
the structure of the outflow inside of r c , which sets the rela¬ 
tive ratio of the outflow velocity to the binary escape veloc- 
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ity, Voo/v esc ~ 0.25 (Fig. 3), the relation between v esc and the 
heat deposition rate, Av/v esc ~ 0.5/p, and the radius r c where 
the heat is deposited. It is reasonable that binary inspiral 
inside a common envelope will also lead to an equatorially- 
concentrated outflow, with a velocity and heating rate near 
its inner edge which are proportional to the instantenous es¬ 
cape speed of the binary (e.g. Taam & Bodenheimer 1989; 
Livio et al. 1990). Naturally, the precise numerical values 
of Voo/vesc, Av/v e sc, r c /a, and even the functional form of 
Lthin(Vesc) will likely be different, and potentially varying in 
time, but the observable features will be described by a sim¬ 
ilar semi-analytic theory. Specifically, the shock heating of 
the ejected material might happen much closer at r c =» a. 
This should not change the outcome qualitatively, because 
a large portion of the energy will be lost to adiabatic expan¬ 
sion before the material can efficiently radiate, as shown in 
our work. 


5 DUST FORMATION 

A necessary condition for dust formation is that the outflow 
temperature decrease below the dust condensation tempera¬ 
ture, r dust . Although this is the only condition to form dust 
within the lowT_f a05_gn93 opacity table (where we esti¬ 
mate r dust ~ 1200 K from the opacity increase; Fig. 2), this 
condition is not sufficient. Following Kochanek (2011b), the 
maximum size of dust grains, a max , which can form at a given 
radius r dust is 


^max 


k'dustW 'll 

167t(z max / r dust)Pbulk r’coCduM 


(35) 


where v dust » lkms -1 is the collision velocity between co¬ 
agulating particles for gas of temperature ~ r dust ~ 10 3 K, 
X g = 0.005 is the mass fraction of condensible material, 
Pbulk ~ 3gem -3 is the bulk density of the grain, z max is the 
vertical thickness of the outflow, and Voo ~ 0.25 v e sc- The con¬ 
dition to form dust is 4a^ lax p du ik/mo > N dust , where niQ = 12m p 
and A dust ~ 7 is some critical number of molecules to start the 
coagulation process. This corresponds to a mass-loss rate of 


M > 10 -9 M 0 yr _1 



(36) 


Equation (36) shows that cool L 2 outflows can in principle 
satisfy the condition for dust grain growth, even at very 
low M. This conclusion is relatively insensitive to the semi¬ 
major axis of the binary. We are thus justified to use the 
lowT_f a05_gn93 opacity tables, which include dust grains. 

We have discovered another effect that can suppress 
dust formation in our simulations. In Figure 17, we show 
the temperature as a function of radius and time for the 
SG case with M = 1CT 5 M 0 yr -1 . The heat generated by the 
spiral shocks is initially efficiently redistributed radially by 
radiative diffusion, because M d i d dominates the energy equa¬ 
tion. This results in the outflow being nearly isothermal for 
r > r c =s 8o (f = 5P in Fig. 17) . As the outflow evolves, however, 
a larger fraction of the power is lost to adiabatic expansion, 
causing the temperature of the isothermal outflow to gradu¬ 
ally decrease with time ( until t = \1P in Fig. 17). When the 
outflows approaches the dust formation temperature, k dra¬ 
matically increases and radiative cooling becomes the domi¬ 
nant term in the energy equation. The outflow temperature 



Figure 17. Evolution of a radial dependence of temperature for 
the binary from Figure 4 with M = 10 5 M 0 yr _1 for five epochs af¬ 
ter the start of the simulation. The shaded and hatched regions 
denote the 5% and 95% of the temperature distributions as a 
function of radius. The inset plot shows time evolution of the lu¬ 
minosity and effective temperature produced by the outflow. For 
clarity, this calculation does not include the irradiation by the 
central binary. 


suddenly decreases as a result, essentially on the diffusion 
timescale (t = 17.8P), and a new equilibrium corresponding 
to a dusty outflow is quickly established (t = 18. 3P). 

The inset plot illustrates this behavior by showing the 
time evolution of the radiative luminosity and effective tem¬ 
perature. The luminosity initially rises as the shocked out¬ 
flow is established, but then asymptotes to a constant value 
as the outflow becomes cooler as described above. When the 
T ~ T dust and dust first forms, the luminosity spikes, releas¬ 
ing the accumulated thermal energy. There is an associated 
decrease in the effective temperature, because cooler areas 
of the outflow now contribute to the luminosity as a result 
of the higher value of k. 

Figure 18 further illustrates this scenario by showing 
the temperature and density evolution of the outflow on a 
background shaded by the optically-thin cooling rate. After 
leaving the L 2 , the gas expands freely along trajectories of 
decreasing temperature and density. Upon encountering the 
previous wrapping of the spiral, it is heated to the tempera¬ 
ture T sd (Eq. [25]). Depending on the binary parameters (M, 
a) and the outflow density, the density and temperature can 
land in the opacity valley (3.2 < logT < 3.6, p < 10~ 9 gem 3 ), 
where the optically-thin cooling rates are several orders of 
magnitude smaller than in the surrounding areas. 

The energy equation (Eq. [8]) is dominated by radiative 
diffusion, which drives the outflow to constant temperature 
at r > r c . As the gas loses energy to adiabatic expansion (and 
to a lesser degree also to radiative cooling), this radially- 
constant temperature decreases with time (Figure 17). The 
rate of adiabatic energy loss for r> r c , where poc r~ 2 , is 


ii 


gas ~ 


■ 2 - — , 
p r 


(37) 


where Voo °c v esc is the asymptotic velocity (assumed to be 
constant for r > r c ). For an ideal gas with P = k^pT/m p , the 
total adiabatic energy loss between r c and the outer radius 


MNRAS 000, 000-000 (2015) 
























16 Pejcha et al. 


lowT fa05 gn93 



2.5 3.0 3.5 4.0 

log T [K] 


Figure 18. Schematic illustration of the origin of the two types 
of outflows (Fig. 10 and 17) and their transition in temperature- 
density space, with the optically-thin cooling rate kctsbT 4 as a 
shaded background. The arrows represent particle densities and 
temperatures from snapshots of our simulations. Initially, the 
gas cools quasi-adiabatically along lines of decreasing temper¬ 
ature and density, until being heated when the encounter the 
shock. If the shock heating is weak, then the gas remains at 
low temperatures and high opacities; dust forms and the result¬ 
ing increased opacity and radiative cooling rate drive the out¬ 
flow to even lower temperatures. Alternatively, if shock heating is 
stronger and places the gas into the p — T region of low radiative 
cooling efficiency, then radiative diffusion dominates and drives 
the outflow to a radially constant temperature. Due to adiabatic 
losses, the temperature of the isothermal outflow also slowly de¬ 
creases as the outer radius r out grows in time. Once the dust for¬ 
mation temperature is reached, the entire outflow abruptly tran¬ 
sitions to the former dust-dominated state. See the text for more 
details. 


of the outflow, r ou t, is thus given by 


J 'Vout 

r c 


M 

^gas — dr — 



m p 



(38) 


The total available power compensating for adiabatic losses 
is temporally constant and equal to Lthin (Eq. 29]), provid¬ 
ing a relationship between the outflow temperature T , r out , 
and the velocity spread of the spiral Av 2 . Once T —> T & ust , ra¬ 
diative cooling suddenly becomes much more efficient, and 
the outflow radiates its thermal energy on the radiative dif¬ 
fusion timescale. This happens once the outflow raches a 
critical size 


1 / r 0 ut\ 1 nipAv~ 4 T sd 

\ r c J 4 A.'s /'diist 3 r (lusl 


(39) 


For higher values of T s h, an exponentially larger outer radius 
is required to reach the dust formation temperature 7)]^. 
For the example shown in Figure 17, T s h ~ 1900 K (t = 5 P in 
Fig. 17) and r dust » 1200 K, which gives r out ~ 8r c . The first 
gas element requires a time t~ 30 P to reach this r ou t, a result 
in broad agreement with Figure 17. Note that because r ou t is 
exponential in T s h, small changes in the binary parameters 
results in dramatic changes in the time of the dust transition 
and the amount of energy released. 

The total amount of energy released when dust forms 


is given by 


3 Fdust M 


-r/aex p 


10 41 ergfT^W 

°l 1200K/ inAs) 


4 T sh \ 
3 7 dust / 

3/2 


r out / r c 
10 


The radiative diffusion timescale across the outflow is 
-2 


tdiff ■ 


it K P 


16s 


r out 

Io7 


5* 


10 4 cm 2 g 1 


P-10. 


(40) 


(41) 


where p - p_iolO~ 10 gem -3 and we assumed typical opac¬ 
ity for the opacity minimum. As the radiative diffusion 
timescale is short in all cases of interest, the duration of 
the luminosity spike is limited by the time to radiate the 
thermal energy, 


u . 

^cool — . ~ 10 s 

M cool 


/ ^dust \ ^ l _^_ 

V 1200K/ \10- 1 cm 2 g- 1 


(42) 


where we used typical opacity for dust grains. Alternatively, 
the finite time required to form dust grains could be the rate 
limiting step in the rapid transition (eq. [36]). The double- 
peaked structure in the time evolution of L in Figure 17 is 
caused by changes in k as T decreases (Fig. 18). 

Note also that the transition between radiation diffu¬ 
sion dominated outflow and dust formation will appear dif¬ 
ferently from different viewing angles. As dust formation 
typically starts at r out and then propagates inward, the bi¬ 
nary will become optically-thick in the orbital plane almost 
immediately, causing the luminosity to drop significantly. 
However, the outflow can remain optically-thin in the ver¬ 
tical direction, where all the radiation escapes, and so after 
dust forms the binary will continue to be bright when viewed 
approximately pole-on. 

There are reasons to believe that the sudden formation 
of dust observed in our simulations might be an artifact of 
the simplifications made in this work. First, radiative diffu¬ 
sion should drive the system to a constant radiative flux (a 
gradient of radiation energy), rather than radially constant 
temperature. However, the L 2 outflow is not thermally sup¬ 
ported and has a finite outer boundary, because the particles 
are being injected into an empty space. When radiative cool¬ 
ing is inefficient, radiative flux cannot leave the system at 
the outer boundary, so the only way to have radially con¬ 
stant flux is with zero flux, i.e. a radially constant tempera¬ 
ture. Also, in physically realizable situations, the L 2 outflow 
does not begin instantaneously, but instead grows on many 
dynamical timescales of the binary, as was potentially ob¬ 
served in V1309 Sco (Pejcha 2014). In such cases, enough 
previous matter might be already present to enable radia¬ 
tion to quickly diffuse out, which will then not present an 
obstacle to dust formation. 

Second, we do not differentiate between the Planck 
and Rosseland opacity means in the expressions for ra¬ 
diation cooling. In the region of low opacity, the Planck 
means might be noticeably higher than the Rosseland means 
(Tomida et al. 2013, Fig. 27). As a result, optically-thin ra¬ 
diative cooling will not be as strongly suppressed and might 
dominate the energy equation, again suppresing the diffusion 
and leading to a rapid dust formation. Finally, the properties 
of the dust such as the opacity depend on the metallicity of 
material, which we assumed to be constant and equal to ap¬ 
proximately solar. Furthermore, the composition of dust in 
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stellar mergers might be very different from the normal dust 
in terms of isotopic ratios (Kaminski et al. 2015). To sum¬ 
marize, our results indicate only the possibility rather than 
a certainty of sudden dust formation. This issue should be 
illuminated with simulations employing more sophisticated 
radiation transport, including more appropriate opacity ta¬ 
bles. 


6 SUMMARY OF OBSERVABLE PROPERTIES 

AND CONNECTION TO RED TRANSIENTS 

We have shown that the luminosity generated by the par¬ 
tial dissipation of the kinetic energy of the outflow from the 
L 2 point is directly proportional to the mass-loss rate and 
the escape velocity of the binary v esc = -\J2GM/a. We have 
also argued that similar observational characteristics will be 
exhibited by a much wider class of models motivated by in¬ 
spiral of a binary star within a common envelope: disk-like 
outflow with asymptotic velocity and energy deposition rate 
at the inner edge proportional to v esc . We now use our results 
to predict the observable characteristics of optical/infrared 
transients produced by binary mass loss. 

The primary characteristics of a transient event are its 
peak luminosity and characteristic duration f 0 t, s (Kasliwal 
2011, 2012). If angular momentum loss through the L 2 point 
dominates the evolution of the merging binary, the timescale 
of the associated transient is set by the response of the en¬ 
velope of the primary star to mass loss and the associated 
angular momentum loss. At later phases, drag forces in the 
non-corotating common envelope will determine the inspiral 
and the energy deposition rates. These processes are chal¬ 
lenging to predict with confidence, thus leading us to instead 
consider a range of timescales as a free parameter and ask 
what is the maximum luminosity produced by a particular 
binary on this timescale. We assume that the binary uni¬ 
formly loses 10% of its mass over the transient duration t a t> s , 
resulting in a the mass-loss rate 

M 

M = 0.1-. (43) 

Gbs 

We also assume that the binary parameters are approx¬ 
imately constant during the transient. In reality, M and 
the binary parameters must evolve during the event, which 
will affect the observed characteristics. Furthermore, the ob¬ 
served properties of the transient are expected to be view¬ 
ing angle dependent, which we discuss later in more detail. 
Our goal here is to show that our results can broadly match 
some of the general characteristics of the red transients, and 
we defer the application our results to specific transients 
to forthcoming work. Thus, we estimate the transient peak 
luminosity as L(f 0 bs)i where L(t) is evaluated using the semi- 
analytic model (Sec. 4). 

In Figure 19, we show the resulting maximum luminos¬ 
ity as a function of the timescale for a grid of binaries bina¬ 
ries with 0.005 < a < 0.5 AU and 1.4 < M < 14 M Q . The lower 
limit of a is set approximately by the size of the core of 
evolved massive stars. For each set of binary parameters we 
investigate f Q b s between 5 and 800 days, which roughly corre¬ 
sponds to the observed timescales of red transients. We show 
only results in cases where the peak transient luminosity is 
at least 10 times higher than the intrinsic luminosity of the 


binary L«. We see that for each binary, the peak luminosity 
initially increases as t Q b s gets longer, because the luminos¬ 
ity is limited by the time it takes the radiation to diffuse 
out of the outflow, which is set by the expansion timescale 
of the outflow (t obs < f S cale)- Eventually, the peak luminosity 
starts decreasing with further increases in r 0 b s , because the 
mass-loss rate corresponding to the given timescale cannot 
produce sufficiently high luminosity (t 0 b s <; f sca i e ). The indi¬ 
vidual lines are relatively well-ordered by the escape velocity 
of the binary with more compact binaries (smaller M/a) pro¬ 
ducing higher luminosities. 

Figure 19 shows that the peak luminosity reachable 
within our framework is — 10 6 Lq. The transients occupy the 
parameter space between dwarf novae, FU Ori outbursts, 
core-collapse supernovae and the bulk of classical novae. 
Note that our results hold even for mechanisms for produc¬ 
ing outflows from binaries other than L 2 mass loss, unless 
the model parameters Av/v e sc an d Voo/v es c differ by more than 
order of magnitude from what we assumed. 

The predictions of our semi-analytic model match the 
luminosities and effective temperatures of the observed class 
of red transients, with the exception of SN2008S, NGC300 
OT, and probably also PTFlOfqs. Reaching luminosities 
necessary to explain these events would require noticeably 
higher values of v eS c, implying a very small semi-major axis. 
This might not be realistic for transients involving massive 
stars such as SN2008S and NGC300 OT, which had progen¬ 
itor masses of 6 -IOMq (Prieto et al. 2008, 2009). Note also 
that the timescales of optimal energy release of each binary 
(maxima of the lines of the semi-analytic model in Fig. 19) 
decrease as v e sc increases, which implies that the predicted 
timescales might be too short if one attempts to reproduce 
the luminosity of SN2008S and NGC300 OT by increasing 
Vesc- Leaving the luminosities and timescales aside, SN2008S 
and NGC300 OT differ in a number of other fundamental 
aspects from the rest of RT, as we detailed in the Introduc¬ 
tion. 

In Figure 20, we show the predicted effective tem¬ 
peratures of our semi-analytic model in comparison with 
those of the known red transients. Our model generally pre¬ 
dicts transients with low temperatures, T e g <; 3000 K, ex¬ 
cept the fastest and brigtest objects, which might reach up 
to ~6000K. Note that our semi-analytic model underpre¬ 
dicts the effective temperature for T e jf > 3000 K, which corre¬ 
sponds to the upper left part of Figure 20. Our T e g roughly 
agree with most red transients, again except for SN2008S 
and NGC300 OT. We also predict that longer transients will 
have lower temperatures, potentially reaching below 1000 K, 
as was observed in OGLE 2002-BLG-360. The effective tem¬ 
perature of CK Vul is unknown, but might have been higher, 
because the object exhibited three faster brightenings over 
~ 700 days. The very low effective temperatures allowed in 
our model are a potential distinguishing feature from pure 
hydrogen recombination models (Ivanova et al. 2013b). 

We predict transient luminosities that decrease rather 
gradually with increasing timescale, which suggests the pos¬ 
sibility ~ 10 5 L o transients with durations of several years 
with r eff < 1000 K. As a result of their low optical luminosi¬ 
ties, such transients may best be discovered by surveys of 
nearby galaxies at near-infrared wavelengths. Indeed, such 
a class of transients may have recently been discovered by 
Spitzer (Kasliwal et al. 2014). 
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Figure 19. Peak luminosities and durations of astronomical transients after Kasliwal (2011, 2012). In addition to objects and classes 
in Kasliwal (2011, 2012), we show dwarf novae (Warner 2003), FU Ori outbursts (Hartmann et al. 1989; Rodriguez et al. 1990), and red 
transients V4332 Sgr (Tylenda et al. 2005a), V1309 Sco (Tylenda et al. 2011), M31 LRN 2015 (Williams et al. 2015), OGLE 2002-BLG- 
360 (Tylenda et al. 2013), and CK Vul (Kato 2003). Colored lines are predictions for mass-losing binaries from our semi-analytic model, 
where each line corresponds to a single set of binary parameters and a range of timescales (left as a free parameter). The color marks 
the escape velocity of the binary, as denoted in the colorbar on the right. 


We expect the spectra of the transients to show 
cool photospheres appropriate for their effective tempera¬ 
tures with emission lines of Balrner series, and potentially 
rotational-vibrational transitions of H 2 and CO in the near- 
infrared. If dust is reprocessing the radiation, we expect low- 
temperature continuum with dust emission features. Some 
features of the infrared spectra might be similar to innermost 
regions of proto-stellar disks due to similarities in densi¬ 
ties and temperatures (e.g. Carmona et al. 2011; Beck et al. 
2012 ). 

In our model the asymptotic velocity of the outflow, as 
well as its heating rate from inter-stream shocks, increases 
with the binary escape velocity. This behavior results nat¬ 
urally from L 2 mass loss (Figs. 3 and 19), but it may ex¬ 
tend more generally to other mechanisms of binary mass 
loss powered by orbital decay. This scaling also results in 
a correlation between the observed luminosity and the ex¬ 
pansion velocity of the outflow. In Figure 21, we show the 
asymptotic outflow velocities of our model in comparison 
to the observed expansion velocities of red transients, as 
typically derived from full-width-half-maximum of Balmer 
lines. The observed transients exhibit rough correlation be¬ 
tween these two quantities, which to our knowledge was not 


demonstrated before. The velocities predicted by our model 
roughly parallel the best-fit line through the data, although 
the model is systematically higher. In some cases, this can be 
attributed to inclination effects, when the disk-like outflow 
is not observed exactly edge-on. In V1309 Sco, the binary 
was observed nearly edge-on (Tylenda et al. 2011), and thus 
part of the differences might be attributed also to uncertain¬ 
ties in our model, differences in the extraction of expansion 
velocity, and correction for the instrumental resolution of the 
spectrograph, which is crucial for low-resolution spectra. 

Ivanova et al. (2013b) proposed that the red transients 
are powered by the hydrogen recombination energy from an 
expanding shell ejected from the merging binary. Analytic 
estimates of the properties of transients produced by dif¬ 
fusion of thermal energy through a recombination front, as 
is commonly applied to Type II-Plateau supernovae (e.g. 
Arnett 1980; Popov 1993; Kasen & Woosley 2009), predict 
that the ejecta velocity should be proportional to the ex¬ 
plosion energy £ e xp and ejecta mass M e j, Voo ~ (E exp /M e j) 1 ^ 2 . 
For a supernova, the explosion energy is set by the explo¬ 
sion mechanism in the inner ~ 1000 km of the star and hence 
there is a great difference between Voo and v e sc- For exam¬ 
ple, assuming that Voo of a supernova is a fraction of v eS c 
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Figure 20. Effective temperatures of the transients pro¬ 
duced by the L 2 mass loss as a function of the peak lu¬ 
minosity and the transient timescale. We also show observed 
red transients labeled by their peak effective temperature in 
Kelvin, including V1309 Sco (Tylenda et al. 2011), V4332 Sgr 
(Tylenda et al. 2005a), PTFlOfqs (Kasliwal et al. 2011b), M85 
OT (Kulkarni et al. 2007), SN2008S (Smith et al. 2009), NGC300 
OT (Berger et al. 2009), V838 Mon (Tylenda 2005), OGLE 2002- 
BLG-360 (Tylenda et al. 2013), M31 RV (Mould et al. 1990), 
and M31 LRN 2015, where we converted the F5I spectral type 
(Kurtenkov et al. 2015) to temperature using observations of 
van Belle et al. (1999). The effective temperature for CK Vul was 
not measured. 


from the surface of a red supergiant would give completely 
misleading results. This illustrates the difference between an 
explosion and a wind-like outflow. For red transients the as¬ 
sumption, that Voo “c v es c gives reasonable results. Both L 2 
mass loss model presented here, and the model proposed by 
Ivanova et al. (2013b), yields a correlation between Voo and 
peak luminosity. 

Kochanek et al. (2014) found that the peak luminos¬ 
ity of the transient is a steep function of the progenitor 
mass, LocM 2 ~ 3 . We would expect more massive progenitors 
to produce more luminous events for both our model and 
the hydrogen recombination model. The exact scaling with 
mass will depend on a combination of ejected mass, ejec¬ 
tion timescale, the diffusion time through the ejecta, and 
the relation between mass and size of the binary. 

We note that, due to the disk-like geometry of the out¬ 
flow, which we have characterized quantitatively for L 2 mass 
loss (Fig. 6), the transient properties may depend on the 
viewing angle of the observer with respect to the orbital 
plane. As M increases, the outflow becomes optically-thick 
first in the orbital plane and, if dust forms, radiation will es¬ 
cape primarly in the perpendicular direction. If the orbital 
plane is viewed face-on, dust acts primarily to decrease T e ff 
and to prolong the timescale for radiation to escape the out¬ 
flow. If dust forms during the transient, potentially through 
the mechanism described in Section 5, then observers view¬ 
ing the event close to the orbital plane will observe fading 
and dust formation. However, face-on observers may observe 
little change in L and T e g, or even a modest brightening 
as rapid cooling releases the accumulated thermal energy 
(Fig. 17). 


Finally, it is clear that the full range of phenomenology 
exhibited by red transients are more complex than can be 
accommodated by our model, and may involve multiple mass 
ejection mechanisms with a geometry more complex than 
the equatorial outflow heated at the inner edge studied here 
(e.g. Thompson et al. 2009; Kochanek 2011a; Nandez et al. 
2014; Kashi & Soker 2015). 


7 DISCUSSION AND CONCLUSIONS 

We have analyzed the radiation hydrodynamics of unbound 
outflows originating from the L 2 point of binary stars, em¬ 
ploying realistic equation of state and opacities. The out¬ 
flowing spiral stream spreads and the spiral arms overlap 
and merge at a collision radius, where the resulting shock 
wave thermalizes ~ 10-20 per cent of the radial kinetic en¬ 
ergy of the outflow (Figs. 4 and 7). Depending on the binary 
properties and mass loss rate, radiative cooling, diffusion, or 
adiabatic expansion dominate the energy evolution of the 
outflow (Fig. 4) with effects on the produced luminosity 
(Fig. 8) and dust formation (Fig. 17). 

We generalize our results to any predominantly equa¬ 
torial outflow originating from a binary star, where the en¬ 
ergy deposition at the inner edge and the asymptotic ve¬ 
locity are both proportional to the binary escape velocity. 
This is motivated by a some scenarios of the common enve¬ 
lope phase. We construct a semi-analytic model of such out¬ 
flows, including luminosity and effective temperature evo¬ 
lution (Sec. 4), and investigate the peak luminosities as a 
function of the timescale (Fig. 19). The produced luminosi¬ 
ties (Fig. 19), effective temperatures (Fig. 20), and ejecta 
velocities (Fig. 21) are able to match red transients such 
as V838 Mon, V1309 Sco, and OGLE 2002-BLG-360, but 
we cannot reasonably match the properties of SN2008S, 
NGC300OT, and PTFlOfqs indicating a different physical 
process for these transients. 

Our results are similar to the models based on hydrogen 
recombination in an expanding shell (Ivanova et al. 2013b; 
Nandez et al. 2014) in the sense that the pseudophotosphere 
in our models is often positioned at the greatest opacity gra¬ 
dient corresponding to hydrogen ionization/recombination 
and that the ejecta velocity is proportional to the binary es¬ 
cape velocity. However, we naturally take into account much 
wider range of densities and temperatures and implicitly in¬ 
clude dust formation. Our model can produce r e ff < 3000 K, 
which would be hard for a purely hydrogen recombination 
model, and we can explain the ~ 850 K effective temper¬ 
ature of OGLE 2002-BLG-360 (Tylenda et al. 2013). Fur¬ 
thermore, we suggest the possibility of bright (L ~ 10 5 Lq), 
cool (r eff < 1500 K), and long (t a b s ~ 1000 days) transients, 
which would radiate predominantly in the near-IR. Such 
transients might have been recently discovered with Spitzer 
(Kasliwal et al. 2014). 

Although we are able to match the general properties 
of the observed class of red transients, the evolution of the 
binary and its orbit is not directly coupled to the outflow in 
our models. As a result, we cannot reproduce the light curve 
of individual transients in detail. For example, V838 Mon 
showed evidence of multiple mass ejections, possibly inter¬ 
acting with each other (Tylenda 2005). Similar complicated 
behavior is observed in many other red transients. Evolving 
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together the hydrodynamics of binary interactions and radi¬ 
ation hydrodynamics of the resulting outflow is challenging 
within a single code due to the large range of timescales 
and physical conditions. Instead, the dynamical results of 
mass ejection from simulations such as those of Nandez et al. 
(2014) could be used as an input for our radiation hydrody¬ 
namics model in future work. 

Explaining portions of the light curves of some red 
transients might be feasible even without detailed model¬ 
ing of the binary hydrodynamics. V1309 Sco exhibited a 
slow and gentle ~ 150 day rise to the maximum, which can 
be reconciled with the preceding orbital period evolution by 
considering the L 2 mass loss, as argued by Pejcha (2014). 
The rate of mass loss from L 2 depends primarily on the 
physics of the stellar surface layers, and hence a simple 
power-law runaway model might be sufficient to describe 
the system and to provide useful constraints on the the¬ 
ory (Paczynski & Sienkiewicz 1972; Webbink 1977; Pejcha 
2014). OGLE 2002-BLG-360 is similar to V1309 Sco in these 
aspects (Tylenda et al. 2013). Investigations of these issues 
with our detailed radiation hydrodynamics models will be 
the subject of future work. 

Another extension of the presented work we are cur¬ 
rently pursuing are the optical signatures of the formation 
of an excretion disk. This should be the outcome of L 2 mass 
loss for binaries with mass ratios q < 0.06 or q > 0.8. For these 
binaries, the outflow remains bound and stalls at a particular 
radius (Shu et al. 1979), forming a mass loss ring that con¬ 
tinues to interacts with the surface of the binary star. Differ¬ 


ent inner boundary conditions will be needed to properly in¬ 
vestigate the structure of the mass-loss ring and its effect on 
the central binary. In a preliminary exploration of the early 
phases with our current inflow boundary condition, we find 
that such systems can be much more luminous and more effi¬ 
cient in extracting angular momentum from the binary than 
in the unbound outflow case. However, a detailed investiga¬ 
tion is beyond this scope of this paper and is relegated to 
future work. Our results might be applicable also to planet- 
star mergers, which previous work has suggested could give 
rise to red-nova like transients (Metzger et al. 2012), and to 
double white dwarf mergers, which might be the common 
outcome of their interaction (Shen 2015). 

Finally, we mention the possibility that the shock heat¬ 
ing of the ejecta might not be stationary at the inner edge of 
the outflow, but the shock wave might dynamically propa¬ 
gate through the outflow. This may be a result of a collision 
of two shells launched with different asymptotic velocity, 
or due to instantenous deposition of a large amount of en¬ 
ergy at the base of the outflow (e.g. Sandquist et al. 1998; 
Morris & Podsiadlowski 2006). In both cases, the propa¬ 
gating shock wave will heat the outflow and the resulting 
electromagnetic emission will depend on the density and 
velocity profile of the outflow, and the energy deposited. 
When the shock reaches the photosphere the radiation be¬ 
hind the shock might be suddenly released in a form of 
a shock breakout. The formalism developed for the shock 
breakout from an optically-thick wind in the context of 
interaction-powered supernovae could be applied to this sit- 
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uation (e.g. Balberg & Loeb 2011; Chevalier & Irwin 2011; 
Ginzburg & Balberg 2012, 2014; Svirski et al. 2012), how¬ 
ever the strong asphericity of the binary outflow will likely 
alter the transient characteristics. 
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APPENDIX A: DEFINITION OF QUANTITIES 
FOR RADIATIVE DIFFUSION, COOLING, AND 
HEATING 

Al Radiative diffusion 


The radiative diffusion coefficient ki for particle i is calcu¬ 
lated as 


ki = 


16cr SB 

pm 


* iTf, 


(Al) 


where Kj is the Rosseland-mean opacity and A, is the flux 
limiter 


2 + R, 

6 + 3R t +Rr' 


(A2) 


The quantity R, is defined as 

„ |V«r,;| 

K-i — •> 

UrjPm 


(A3) 


where we calculate the gradient of the radiation energy den¬ 
sity u t j = 4 cr SB Tpc as 


VM r J - 



Ur XWh. 


(A4) 


The calculation of radiative diffusion is very fast, because it 
depends only on local quantities. 


A2 Radiative cooling and irradiation by the 
central star 

We estimate the column density T, z j and optical depth t z j 
by integrating p and xp, respectively, along a ray going per¬ 
pendicularly out of the orbital plane and originating at a 
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Figure Al. Heating rate due to the irradiation from the central star for a calculation with M = 10“ 3 M©yr -1 , a = 0.03 AU, M 2 = 1.5 M©, 
q = 0.15, and N = 3000. Left panel shows the radial dependence of E r j for 20 randomly chosen particles. Right panel shows the heating rate 
Wheat as a function of position. Each point is a single particle and the color represents Wheat as indicated by the color bar. The surface of 
the binary (red closed shape near the center) is replaced by a point source at the barycenter (red star). Following the material as being 
ejected from L 2 , the highest Wheat occurs on the inner rim of the stream and the outer rim is partially self-obscured. The heating rate 
drops substantially after the spiral starts a second winding and then Wheat only gradually increases as the interior material expands and 
becomes more transparent. 


particle j. Specifically, we use the SPH estimator to evalu¬ 
ate p and Kp along this ray as 

s ;j= ffldz- V mi f W(\rj-ri\,hj)dz, (A5) 

J i,Zi/zj> 1 

T Z J= f/cpdz~ 2 m,K, f W(lrj-ri\,hj)dz, (A6) 

J i,Zi/Zj> 1 J ~°° 

where the sums over the particles includes all particles i, 
which are located farther out of the orbital plane than the 
particle j, Zilzj > 1. We do not take into account the con¬ 
tributions to the optical depth from particles with Zi/zj < 1, 
which lie in the smoothing sphere of j. We assume that all 
particles with nlzj > 1 contribute to the column density and 
optical depth fully, so that we can integrate over the full 
smoothing volume. The integral over the smoothing kernel 
can be expressed analytically, but for the sake brevity we do 
not reproduce the full expression, and in our calculation, we 
precompute the values in a lookup table. For the purposes 
of finding the nearest neighbors of particles along the ray 
perpendicular to the orbital plane, we project the particles 
on the orbital plane and construct a two-dimensional tree. 
This brings a relatively small overhead, because the particles 
do not get very far from the disk and the hence the number 
of nearest neighbors is not significantly higher than in three 
dimensions for the same smoothing length. 

To estimate the effects of irradiation of the central star 
(Eq. [13]), we calculate the radial column density and 
optical depth T r j as in Equations (A5-A6), with the ex¬ 
ception that we are interested in a radial ray. To find the 
particles along the radial ray, we project all particles on a 
unit sphere, construct a three-dimensional tree on this pro¬ 
jection, and find all nearest neighbors within the smoothing 
length hj!rj, effectively obtaining particles in a radial cone. 
Although the actual construction of the three-dimensional 
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Figure Bl. Luminosity and effective temperature as a function 
of the stream spreading parameter e for the binary parameters 
from Figure 4, and N = 3000/P. 


tree is relatively fast, the number of particles in a given 
cone can be relatively large, because the smoothing length 
typically increases linearly with the radial distance. As a 
result, calculating Wheat,;' is relatively expensive and signifi¬ 
cantly slows down the calculation. To mitigate this, we in¬ 
clude only particles with rj < 15fl. Outside of this radius, the 
equilibrium temperature of the star is low and the radia¬ 
tive cooling dominates. In Figure Al, we show the typical 
behavior of E r j and Wheat- 
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Figure B2. Radiative cooling luminosity (left panel) and effective temperature (right panel) as a function of elapsed time for a simulation 
with a = 0.03 AU, M = 0.43 M Q , q = 0.13. and M = 10 3 /V/ l0 yr 1 . The lines are for different rates of particle ejection N = 10 3 /P (black), 3x 10 3 /P 
(red), 10 4 /P (blue), and 3x10 4 /P (green). 


APPENDIX B: DEPENDENCE ON 
TECHNICAL PARAMETERS 

In addition to physical parameters, it is worth investigat¬ 
ing how our results depend on the technical parameters 
of the calculation. Although the stream-spreading param¬ 
eter s should be approximately constant for different binary 
classes (Eq. [20]), the specific value is not known without 
detailed hydrodynamical study of the binary surface at the 
L 2 point. In Figure Bl, we show the dependence of the lu¬ 
minosity and effective temperature on e. We see that for 
e > 0.02 both L and T e ff increase, because wider stream im¬ 
plies greater velocity spread and higher shock temperatures 
(Sec. 4). For e < 0.02 both L and T e ff stay approximately con¬ 
stant. The reason is that for too small e the particles will 
be injected in a much smaller region than what would corre¬ 
spond to their temperature, the thermal pressure forces will 
become important and the stream will expand. Figure Bl 
also shows that our fiducial e = 0.05 perhaps slightly overes¬ 
timates the appropriate value, but the differences are much 
smaller than other uncertainties in the input physics and 
calculation methods. 

In Figure B2, we show the dependence of the luminosity 
and effective temperature on the rate of particle injection N, 
which corresponds to the simulation resolution. The scatter 
of the lines goes down as N increases, especially during the 
transient epoch when the spiral arms start merging for the 
first time. Runs with higher N also converge to the asymp¬ 
totic value faster, but the asymptotic luminosities and effec¬ 
tive temperatures are independent of resolution. We run all 
our calculations for at least 15P, but typically much longer 
until the steady-state is achieved. 
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